{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "03ee1f11-7a60-4761-aaf7-673ea04e7dca",
   "metadata": {},
   "source": [
    "# Polygon Tutorial\n",
    "\n",
    "`h3-py` can convert between sets of cells and GeoJSON-like polygon and multipolygon shapes.\n",
    "\n",
    "We use the abstract base class `H3Shape` and its concrete child classes `LatLngPoly` and `LatLngMultiPoly`\n",
    "to represent these shapes. \n",
    "Any references or function names that use \"H3Shape\" will apply to both `LatLngPoly` and `LatLngMultiPoly` objects.\n",
    "\n",
    "`h3-py` is also compatible with Python objects that implement `__geo_interface__` (https://gist.github.com/sgillies/2217756),\n",
    "making it easy to interface with other Python geospatial libraries.\n",
    "We'll refer to any such object as a \"geo\" or \"geo object\".\n",
    "\n",
    "`LatLngPoly` and `LatLngMultiPoly` both implement `__geo_interface__` (making them geo objects themselves),\n",
    "and can be created from other geo objects, like Shapely's `Polygon` or `MultiPolygon` objects that occur\n",
    "when using `geopandas`.\n",
    "\n",
    "Below, we'll explain the `h3-py` API for working with these shapes, and how the library interacts\n",
    "with other Python geospatial libraries.\n",
    "\n",
    "To start, we'll import relevant libraries, and define some plotting helper functions to visualize the shapes we're dealing with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f5454f0-f65e-491c-b534-f1ef08320263",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h3\n",
    "\n",
    "import geopandas\n",
    "import geodatasets\n",
    "import contextily as cx\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def plot_df(df, column=None, ax=None):\n",
    "    \"Plot based on the `geometry` column of a GeoPandas dataframe\"\n",
    "    df = df.copy()\n",
    "    df = df.to_crs(epsg=3857)  # web mercator\n",
    "\n",
    "    if ax is None:\n",
    "        _, ax = plt.subplots(figsize=(8,8))\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "    df.plot(\n",
    "        ax=ax,\n",
    "        alpha=0.5, edgecolor='k',\n",
    "        column=column, categorical=True,\n",
    "        legend=True, legend_kwds={'loc': 'upper left'},\n",
    "    )\n",
    "    cx.add_basemap(ax, crs=df.crs, source=cx.providers.CartoDB.Positron)\n",
    "\n",
    "\n",
    "def plot_shape(shape, ax=None):\n",
    "    df = geopandas.GeoDataFrame({'geometry': [shape]}, crs='EPSG:4326')\n",
    "    plot_df(df, ax=ax)\n",
    "\n",
    "\n",
    "def plot_cells(cells, ax=None):\n",
    "    shape = h3.cells_to_h3shape(cells)\n",
    "    plot_shape(shape, ax=ax)\n",
    "\n",
    "\n",
    "def plot_shape_and_cells(shape, res=9):\n",
    "    fig, axs = plt.subplots(1,2, figsize=(10,5), sharex=True, sharey=True)\n",
    "    plot_shape(shape, ax=axs[0])\n",
    "    plot_cells(h3.h3shape_to_cells(shape, res), ax=axs[1])\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a41e2b-0c9e-45f5-8005-51af2ef8af32",
   "metadata": {},
   "source": [
    "# LatLngPoly\n",
    "\n",
    "We can create a simple `LatLngPoly` object by providing a list of the **latitude/longitude pairs** that describe its exterior.\n",
    "Optionally, holes can be added to the polygon by appending additional lat/lng lists do describe them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e9e3e2a-a4cf-43df-a498-ee7d15a8cc4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "outer = [\n",
    "    (37.804, -122.412),\n",
    "    (37.778, -122.507),\n",
    "    (37.733, -122.501)\n",
    "]\n",
    "\n",
    "poly = h3.LatLngPoly(outer)\n",
    "print(poly)\n",
    "plot_shape(poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6872b341-afc9-4405-b27a-b5f2b8a19847",
   "metadata": {},
   "outputs": [],
   "source": [
    "hole1 = [\n",
    "    (37.782, -122.449),\n",
    "    (37.779, -122.465),\n",
    "    (37.788, -122.454),\n",
    "]\n",
    "\n",
    "poly = h3.LatLngPoly(outer, hole1)\n",
    "print(poly)\n",
    "plot_shape(poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d7ae0fb-fa96-46fe-8d06-7778575642dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "hole2 = [\n",
    "    (37.771, -122.484),\n",
    "    (37.761, -122.481),\n",
    "    (37.758, -122.494),\n",
    "    (37.769, -122.496),\n",
    "]\n",
    "\n",
    "poly = h3.LatLngPoly(outer, hole1, hole2)\n",
    "print(poly)\n",
    "plot_shape(poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62317cfd-8a0e-4c8a-b899-131acdec4305",
   "metadata": {},
   "source": [
    "## String representation and attributes\n",
    "\n",
    "The `LatLngPoly` string representation given by its `__repr__` shows the number of vertices in the outer loop of the polygon, followed\n",
    "by the number of vertices in each hole.\n",
    "\n",
    "A representation like\n",
    "\n",
    "```\n",
    "<LatLngPoly: [3/(3, 4)]>\n",
    "```\n",
    "\n",
    "denotes a polygon whose outer boundary consists of 3 vertices and which has 2 holes,\n",
    "with 3 and 4 vertices respectively.\n",
    "\n",
    "We can access the coordinates describing the polygon through the attributes:\n",
    "\n",
    "- `LatLngPoly.outer` gives the list of lat/lng points making up the outer loop of the polygon.\n",
    "- `LatLngPoly.holes` gives each of the lists of lat/lng points making up the holes of the polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97c8aa44-12fb-4c47-9254-d931522575ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = h3.LatLngPoly(outer, hole1, hole2)\n",
    "poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "150ae04c-bae2-4b7e-88d0-f5fcf0835e7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly.outer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84061e68-60f3-4df4-a279-a0630af7dcc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly.holes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae175854-6770-44aa-b796-8d9390f7c856",
   "metadata": {},
   "source": [
    "## `__geo_interface__`\n",
    "\n",
    "`LatLngPoly.__geo_interface__` gives a GeoJSON representation of the polygon as described in https://gist.github.com/sgillies/2217756\n",
    "\n",
    "**Note the differences in this representation**: Points are given in **lng/lat** order (but the `LatLngPoly` constructor expects **lat/lng** order), and the last vertex repeats the first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb79c742-27aa-481c-971b-28c100109efe",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = poly.__geo_interface__\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc1b0ae8-bcce-4e5f-ad88-57e3797518f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(poly.outer)\n",
    "print(poly.holes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f567fae-9f4e-472f-afcd-bd6f7bb32cba",
   "metadata": {},
   "source": [
    "We can create an `LatLngPoly` object from a GeoJSON-like dictionary or an object that implements `__geo_interface__` using `h3.geo_to_h3shape()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42d7d580-2bba-4cd3-8d13-2434323e4081",
   "metadata": {},
   "outputs": [],
   "source": [
    "h3.geo_to_h3shape(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70b6fa3-d85d-420e-a295-a75a32bf3dad",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MockGeo:\n",
    "    def __init__(self, d):\n",
    "        self.d = d\n",
    "\n",
    "    @property\n",
    "    def __geo_interface__(self):\n",
    "        return self.d\n",
    "\n",
    "\n",
    "geo = MockGeo(d)\n",
    "h3.geo_to_h3shape(geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ab82119-9d9e-465d-94b5-61083e9d9df3",
   "metadata": {},
   "source": [
    "Also note that `LatLngPoly.__geo_interface__` is equivalent to calling `h3.h3shape_to_geo()` on an `LatLngPoly` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd9bd77d-ed1e-46bb-991a-f0cfc10f634b",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly.__geo_interface__ == h3.h3shape_to_geo(poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55c31b56-c65d-4e06-bf60-edf96bb57e93",
   "metadata": {},
   "source": [
    "## Polygon to cells\n",
    "\n",
    "We can get all the H3 cells whose centroids fall within an `LatLngPoly` by using `h3.h3shape_to_cells()` and specifying the resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c310200-30f2-4b24-a967-d705ffc05bd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "h3.h3shape_to_cells(poly, res=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec0d8ac7-ba97-46c4-8740-3a12a00e318d",
   "metadata": {},
   "source": [
    "We'll use a helper function to show a few different resolutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac7cd51-1c7e-4509-94b4-92e59fd22c2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_shape_and_cells(poly, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0330d60-0adc-46e7-9689-7be1ba957dfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_shape_and_cells(poly, 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fd40b0a-34fd-4ac7-9c3a-c367e5041cd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_shape_and_cells(poly, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faa6a726-2adc-40cc-9e1e-5d4ee66a08ea",
   "metadata": {},
   "source": [
    "## H3 Polygons don't need to follow the right-hand rule\n",
    "\n",
    "`LatLngPoly` objects do not need to follow the \"right-hand rule\", unlike GeoJSON Polygons. \n",
    "The right-hand rule requires that vertices in outer loops are listed in counterclockwise\n",
    "order and holes are listed in clockwise order.\n",
    "`h3-py` accepts loops in any order and will usually interpret them as the user intended, for example,\n",
    "converting to sets of cells. However, `h3-py` won't re-order your loops to\n",
    "conform to the right-hand rule, so be careful if you're using `__geo_interface__` to plot them.\n",
    "\n",
    "Obeying the right-hand rule is only a concern when creating `LatLngPoly` objects from external input; `LatLngPoly` or `LatLngMultiPoly`\n",
    "objects created through `h3.cells_to_h3shape()` **will respect the right-hand rule**.\n",
    "\n",
    "For example, if we reverse the order of one of the holes in our example polygon above,\n",
    "the hole won't be rendered correctly, but the conversion to cells will remain unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34ef8624-c852-419d-9e98-e9ec70ec6bc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Respects right-hand rule\n",
    "poly = h3.LatLngPoly(outer, hole1, hole2)\n",
    "plot_shape_and_cells(poly, res=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e52ae550-e34d-406b-9c82-b283fd34a170",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Does not respect right-hand-rule; second hole is reversed\n",
    "# Conversion to cells still works, tho!\n",
    "poly = h3.LatLngPoly(outer, hole1[::-1], hole2)\n",
    "plot_shape_and_cells(poly, res=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1940071-eb32-44bc-b07c-2a2eef14b21e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Does not respect right-hand-rule; outer loop and second hole are both reversed\n",
    "# Conversion to cells still works, tho!\n",
    "poly = h3.LatLngPoly(outer[::-1], hole1[::-1], hole2)\n",
    "plot_shape_and_cells(poly, res=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0b13d4a-a376-4a8f-b314-898d9cf94c59",
   "metadata": {},
   "source": [
    "# LatLngMultiPoly\n",
    "\n",
    "An `LatLngMultiPoly` can be created from `LatLngPoly` objects. The string representation of the `LatLngMultiPoly`\n",
    "gives the number of vertices in the outer loop of each `LatLngPoly`, along with the number of vertices\n",
    "in each hole (if there are any).\n",
    "\n",
    "For example `<LatLngMultiPoly: [3], [4/(5,)]>` represents an `LatLngMultiPoly` consisting of two `LatLngPoly` polygons:\n",
    "\n",
    "- the first polygon has 3 outer vertices and no holes\n",
    "- the second polygon has 4 outer vertices and 1 hole with 5 vertices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b0b4b2-2cfb-461d-bd69-511c87e208a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly1 = h3.LatLngPoly([(37.804, -122.412), (37.778, -122.507), (37.733, -122.501)])\n",
    "poly2 = h3.LatLngPoly(\n",
    "    [(37.803, -122.408), (37.736, -122.491), (37.738, -122.380), (37.787, -122.39)],\n",
    "    [(37.760, -122.441), (37.772, -122.427), (37.773, -122.404), (37.758, -122.401), (37.745, -122.428)]\n",
    ")\n",
    "mpoly = h3.LatLngMultiPoly(poly1, poly2)\n",
    "\n",
    "print(poly1)\n",
    "print(poly2)\n",
    "print(mpoly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2faf991-9d3b-4e0a-8fd4-8c31e19fd07b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_shape(mpoly)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "850912ae-85f1-470c-aeea-1140f549fc95",
   "metadata": {},
   "source": [
    "## MultiPolygon to cells\n",
    "\n",
    "`h3.h3shape_to_cells()` works on both `LatLngMultiPoly` and `LatLngPoly` objects (both are subclasses of `H3Shape`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "944564c0-8198-49ae-81de-cf980dc4b348",
   "metadata": {},
   "outputs": [],
   "source": [
    "cells = h3.h3shape_to_cells(mpoly, res=9)\n",
    "plot_cells(cells)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3712d99-cbe2-4f42-80e3-0acf3fcb8e25",
   "metadata": {},
   "source": [
    "## LatLngMultiPoly affordances\n",
    "\n",
    "- Calling `len()` on an `LatLngMultiPoly` gives the number of polygons\n",
    "- You can iterate through a `LatLngMultiPoly`, with the elements being the underlying `LatLngPoly`s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "217d6e5a-bb8b-4923-b311-0f84ebf48872",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(mpoly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03f6cd31-b18f-4fe7-b259-78741f12c95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in mpoly:\n",
    "    print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eca4cec-2881-4ed1-a65c-3aea38ac9c63",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(mpoly)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fd0a70c-d0cf-487c-b492-6b459dc13a30",
   "metadata": {},
   "source": [
    "## `__geo_interface__`\n",
    "\n",
    "`LatLngMultiPoly` implements `__geo_interface__`, and `LatLngMultiPoly` objects can also be created through `h3.geo_to_h3shape()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2731474c-b3f8-4438-8932-62af21f8d6cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = mpoly.__geo_interface__\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5cc947e-0672-492d-a4bf-f303afbb5cfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = MockGeo(d)\n",
    "h3.geo_to_h3shape(geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b71ee2-deed-4336-8b2c-2ab740555188",
   "metadata": {},
   "source": [
    "# Cells to LatLngPoly or LatLngMultiPoly\n",
    "\n",
    "If you have a set of H3 cells that you would like to visualize, you may want to convert them to `LatLngPoly` or `LatLngMultiPoly` objects using `h3.cells_to_h3shape()` and then use `__geo_interface__` to get their GeoJSON representation. Or you could\n",
    "use `h3.cells_to_geo()` to get the GeoJSON dictionary directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2305120-adf4-4cee-9e2a-58e793e55cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = h3.latlng_to_cell(37.804, -122.412, res=9)\n",
    "cells = h3.grid_ring(h, 2)\n",
    "cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e242238-af2c-45a3-acd0-8daf3773a2bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cells(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9d124ff-4e49-4e34-a6be-4a49bc00bcd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "h3.cells_to_h3shape(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3ca5d18-a78a-4bc6-80d9-d98099da924d",
   "metadata": {},
   "outputs": [],
   "source": [
    "str(h3.cells_to_geo(cells))[:200]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ae6d480-4807-4ca1-ac20-4095a20ded31",
   "metadata": {},
   "source": [
    "# API Summary\n",
    "\n",
    "```\n",
    "\"Geo object\" <-> H3Shape <-> H3 cells\n",
    "```\n",
    "\n",
    "\"Geo objects\" are any Python object that implements `__geo_interface__`. This is a widely used\n",
    "standard in geospatial Python and is used by libraries like `geopandas` and Shapely.\n",
    "\n",
    "`H3Shape` is an abstract class implemented by `LatLngPoly` and `LatLngMultiPoly`, each of which\n",
    "implement `__geo_interface__` and can be created from external \"Geo objects\":\n",
    "\n",
    "- `geo_to_h3shape()`\n",
    "- `h3shape_to_geo()`\n",
    "\n",
    "`H3Shape` objects can be converted to and from sets of H3 cells:\n",
    "\n",
    "- `h3shape_to_cells()`\n",
    "- `cells_to_h3shape()`\n",
    "\n",
    "For convience, we provide functions that hide the intermediate `H3Shape` representation:\n",
    "\n",
    "- `geo_to_cells()`\n",
    "- `cells_to_geo()`\n",
    "\n",
    "## Equivalance notes\n",
    "\n",
    "Note that if an object `a` is an `H3Shape`\n",
    "then the following two expressions are equivalent:\n",
    "\n",
    "- `h3shape_to_geo(a)`\n",
    "- `a.__geo_interface__`\n",
    "\n",
    "Also note that if `a` is an `H3Shape`, then `a` \"passes through\" `geo_to_h3shape()`:\n",
    "\n",
    "```python\n",
    "geo_to_h3shape(a) == a\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5838f6c8-aacb-4813-be7c-756c87945476",
   "metadata": {},
   "source": [
    "# Interfacing with GeoPandas and other libraries\n",
    "\n",
    "The `__geo_interface__` compatibility allows `h3-py` to work with `geopandas` and other geospatial libraries easily.\n",
    "\n",
    "To demonstrate, we start with a GeoPandas `GeoDataFrame` describing the five New York City boroughs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a872e2a1-dd51-4262-bc98-800500121bc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = geopandas.read_file(geodatasets.get_path('nybb'))\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbca69dc-1722-4be4-a1bb-329d7e77ec83",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df(df, column='BoroName')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a060b933-8f32-402b-ac61-7cd04dcf70df",
   "metadata": {},
   "source": [
    "## Use a compatible CRS before applying H3 functions\n",
    "\n",
    "The function `h3.geo_to_cells(geo, res)` takes some \"geo object\" that implements `__geo_interface__` (https://gist.github.com/sgillies/2217756) — like a Shapely Polygon or MultiPolygon — and converts it to the set of cells whose centroids are contained in `geo`.\n",
    "\n",
    "**Common Pitfall**: Be careful about what CRS you are using. `h3-py` expects coordinates as latitude-longitude pairs. The CRS of the current dataframe (EPSG:2263) gives coordinates in feet! You should first convert the data to something compatible. A common choice is **EPSG:4326/WGS84**. You'll get incorrect results if you apply `h3.geo_to_cells` without converting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83fced81-89f1-4795-84b8-c852c8c3f938",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original CRS\n",
    "df.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfaf3b8a-2e53-4adf-87a5-ef9fe616df32",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.geometry[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00e1026a-0a0c-448f-9223-8ce6d4d5a9a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Converting to EPSG:4326/WGS84\n",
    "df = df.to_crs(epsg=4326)\n",
    "df.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81071108-3c13-4e10-8608-892ed3742054",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the `geometry` column now has coordinates in degrees\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01ca62b0-61f3-4204-bb0f-a14c9859cd4d",
   "metadata": {},
   "source": [
    "## Converting a geo to H3 cells\n",
    "\n",
    "First, we select one of the boroughs and get the Shapely `MultiPolygon` describing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0be200e1-2536-4430-90ef-38dc72a75f9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = df.geometry[0]\n",
    "type(geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a3c3c9c-b547-4a69-9030-fb88f486c4aa",
   "metadata": {},
   "source": [
    "Note that Shapely `MultiPolygon` objects implement `__geo_interface__`, so we can use `h3.geo_to_cells()` to get the associated set of H3 cells at various resolutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1a12743-b5e2-4cd1-a7fd-848e920fa928",
   "metadata": {},
   "outputs": [],
   "source": [
    "h3.geo_to_cells(geo, res=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27d24ca5-123f-4462-b37d-bcba6acc4f64",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cells(h3.geo_to_cells(geo, res=9))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00d088e-47b9-4bca-a8b3-d22dab747be3",
   "metadata": {},
   "source": [
    "## Converting all geos in a GeoDataFrame to cells\n",
    "\n",
    "We can apply `geo_to_cells()` to each of the Shapely geometries in the `df.geometry` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9901cc3-cd8d-4543-a8fb-98178010700e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_column = df.geometry.apply(lambda x: h3.geo_to_cells(x, res=8))\n",
    "cell_column"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502efd48-2c0b-4860-8796-7ea2e4e88687",
   "metadata": {},
   "source": [
    "## Converting cells to \"geo objects\"\n",
    "\n",
    "If we assign `df.geometry = cell_column` we'll get an error because the `geometry` column of a `geopandas.GeoDataFrame` must contain valid geometry objects.\n",
    "We can obtain compatible objects by converting the cells to `H3Shape` by applying `h3.cells_to_h3shape()`.\n",
    "\n",
    "(Note that, unfortunately, Pandas has some logic to identify iterable members of a series and then renders a tuple of the elements, rather than our preferred `LatLngMultiPoly.__repr__` representation.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5cb00c-9796-4356-8e2a-8dc50c86103c",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_column = cell_column.apply(h3.cells_to_h3shape)\n",
    "shape_column"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1965d53f-2bba-4840-a627-1bbf55405be3",
   "metadata": {},
   "source": [
    "Note that the column now consists of `LatLngPoly` and `LatLngMultiPoly` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b3eedce-2570-4708-992e-e6d739eecc2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_column[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "014bc63d-5cfa-4a65-95e6-390ddddd1703",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_column[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae6afb2-8267-492a-836e-cf371df91e89",
   "metadata": {},
   "source": [
    "Now, if we assign `df.geometry = shape_column`, our `H3Shape` objects will automatically be converted to Shapely Polygon and MultiPolygon objects via `__geo_interface__`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5cd4cf-226b-405f-9c36-7f1ccc09c342",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.geometry = shape_column\n",
    "df.geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2608ce79-f2c5-4dbc-bf22-a28255343aea",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(df.geometry[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d568b20e-88cb-4269-99ce-dc96a73d93d1",
   "metadata": {},
   "source": [
    "## Visualizing H3 cells\n",
    "\n",
    "We took some Shapely geometry objects, converted them to H3 cells, and converted those back to Shapely geometry objects in a `geopandas.GeoDataFrame`. \n",
    "\n",
    "We can visualize the results with our helper function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ab3ec1-49c2-482d-8b84-154dd8196eda",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df(df, column='BoroName')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
